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ABSTRACT 

We report the results of a series of AMR radiation-hydrodynamic simulations 
of the collapse of massive star forming clouds using the ORION code. These sim- 
ulations are the first to include the feedback effects protostellar outflows, as well 
as protostellar radiative heating and radiation pressure exerted on the infalling, 
dusty gas. We find that outflows evacuate polar cavities of reduced optical depth 
through the ambient core. These enhance the radiative flux in the poleward di- 
rection so that it is 1.7 to 15 times larger than that in the midplane. As a result 
the radiative heating and outward radiation force exerted on the protostellar disk 
and infalling cloud gas in the equatorial direction are greatly diminished. This 
simultaneously reduces the Eddington radiation pressure barrier to high-mass 
star formation and increases the minimum threshold surface density for radia- 
tive heating to suppress fragmentation compared to models that do not include 
outflows. The strength of both these effects depends on the initial core surface 
density. Lower surface density cores have longer free-fall times and thus massive 
stars formed within them undergo more Kelvin contraction as the core collapses, 
leading to more powerful outflows. Furthermore, in lower surface density clouds 
the ratio of the time required for the outflow to break out of the core to the core 
free-fall time is smaller, so that these clouds are consequently influenced by out- 
flows at earlier stages of collapse. As a result, outflow effects are strongest in low 
surface density cores and weakest in high surface density one. We also find that 
radiation focusing in the direction of outflow cavities is sufficient to prevent the 
formation of radiation pressure-supported circumstellar gas bubbles, in contrast 
to models which neglect protostellar outflow feedback. 
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Introduction 



Stars of all masses undergo energetic, bipolar mass loss during their formation (Shepherd 
2003] |Richerlt~aI1[2000| . These outflows feed energy back into large-scale turbulent motions 
that support clouds against collapse, may play a role in dispersing some localized regions 



entirely ( 


Norman & Silk 


1980 


McKee 


1989 


Nakamura & Li 


2007 


Carroll et al. 


2010 


Arce 


et al. 


2010 


) and regulate the final mass of the central star ( 


Matzner & McKee 


2000 


Arce 



& Sargent 2006 Wang et al. 2010). Massive stars likely provide the dominant source of 
radiation feedback in the evolution of their parent molecular clouds and any subsequent 
star formation therein. Although much progress has been made both observationally and 
theoretically, a comprehensive picture of massive star formation and the role of feedback from 



massive stars in mediating the star formation process remains to be elucidated (Zinnecker 



& Yorke 2007). 



Observational evidence supports a picture where accretion and outflow ejection processes 
at work in the formation of high-mass stars proceed as a "scaled-up" version of their low- 
mass, solar-type counterparts. Interferometric molecular line measurements have detected 
quiescent compact cores within dense, infrared dark clouds of mass ~ 100M Q flSwift|2009) as 
likely candidates to the onset of high-mass star formation. Observational surveys have estab- 



lished a correlation between molecular outflow mass-loss and source luminosity (Shepherd 



& Churchwell 



1996 



Richer et al.||2000 ) and between circumstellar mass and luminosity from 



0.1 to 1O 5 L (Saraceno et al. 1996| Chandler fc Richer||2000 ). Several authors have detected 



molecular outflows from massive 


protostars with collimation factors of 2 to 10 ( 


Zhang et 


al. 


2001 


Beuther et al. 


2002alc 




2003, 


2004; 


Qiu et al. 


2007 


Lopez-Sepulcre et al. 


2011 


)• 



similar to that of low- mass stars (Bachiller 1996). Radio thermal continuum emission jets 



commonly associated with low-mass protostars (Rodriguez 1997), have also been identified 
near protostellar sources as luminous as ~ 1O 5 L (Torrefies et al. 1997 Curiel et al. 2006). 
Detection of synchrotron emission arising from the jet in one massive young stellar object 
gives support to the notion of a common magnetic driving mechanism to protostellar outflow 



from stars of all masses ( Carrasco-Gonzalez et al. 2010). Because high- mass accretion disks 



are deeply embedded in dusty envelopes, they are particularly difficult to observe directly. 



A few such detections have, however, been made by maser emission sources (Hutawarakorn 



et al. 2002), in high-resolution submillimeter dust emission (Patel et al. 2005) and in near- 



infrared observations where winds from nearby sources have cleared dust from the line of 
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sight (Niirnberger et al. 2007). These observations suggest that massive stars form through 



disk accretion in direct analogy to the formation of low-mass stars. 

Several key theoretical aspects distinguish high and low-mass star formation despite 
the broad similarity of the observed outflow and ejection phenomena. Massive stars are 
shorter-lived and produce more sources of energetic feedback into their environment than 
their low-mass counterparts. O stars radiate their gravitational binding energy and reach 
the main sequence on Kelvin-Helmholtz timescales of < 10 4 yr whereas solar type stars 
require > 10 7 yr. Stars with masses > 10M Q therefore begin nuclear burning while they are 



still embedded within and accreting from the circumstellar envelope (Shu et al. 1987). The 



resultant spherically-averaged radiation pressure on dust grains in the infalling gas exceeds 



the gravitational pull from the central star (Larson & Starrfield 1971). Massive stars can 



therefore only form by accretion if some mechanism is in place to focus the outward radiative 
flux away from the infalling envelope. A variety of focusing mechanism have been suggested 
QNakano||1989~||Jijina fc Adamsfl9~9~o] |Yorke fc Sonnhalter|[2002| |Krumholz et al.|2009t|McKee 



& Ostriker 2007 Kuiper et al. 2010), including the one of central interest for this paper, 



beaming of radiation in the cavities produced by protostellar outflows. (Krumholz et al. 



2005). Once the embedded stars reach the main sequence, ionizing photons generate HII 



regions, strongly affecting the physical structure and chemistry of their environment. Recent 
observation suggests the existence of stars as massive as 3OOM dCrowther et al.|2010[ ), and it 
remains unclear if such large mass can be reached by accretion alone in spite of these strong 
feedback effects. Massive stars appear predominantly in denser clusters than low-mass stars, 



and massive stars more frequently occur in binary and small-multiple systems (Preibisch et 



aI7||2001 Shatsky fc Tokovinin||2002 Lada||2006 ). Furthermore, recent theoretical (Krumholz 



& McKee 2008), numerical (Krumholz et al. 2010) and observational ( Lopez- Sepulcre et 



al. 2010) evidence indicate a minimum prestellar core surface density for high-mass star 



formation, giving rise to a specific environmental dependence that distinguishes the case of 
massive star formation. 

In this paper we present a series of AMR radiation-hydrodynamic simulations of the col- 



lapse of massive star forming clouds using the ORION code (Truelove 1997 Truelove et al. 



1998 Klein 1999 Krumholz et al. 2007b). These simulations are the first to simultaneously 



include radiation and protostellar outflow feedback, and to study their interaction. This 



work is complementary to that of Peters et al. (2010), which included the effect of photoion 



ization but not of radiation pressure or outflows. To probe the environmental dependence 
for massive star formation, we examine the effect of outflows in star forming cores at several 
surface densities representative of typical massive star forming regions in the Milky Way to 
regions characteristic of extragalactic super star clusters. To isolate the effect of outflow 
feedback alone, we include one model where outflows have been turned off in an otherwise 



-4- 



identical cloud. In ^2] we describe the simulation methodology and input parameters, in £|3] 
the numerical results are presented and discussed and in §|4] we summarize the conclusions 
that can be drawn from the models. 



2. Simulation Setup 



2.1. Initial Conditions 



Our simulations are initialized to the state of a prestellar core of mass M, each with a 
power law density profile given by 

p(r) oc r- kp , (1) 

with k p = 3/2, consistent with models ( McKee fc Tan|[2002 , 2003) and observation (Beuther 
et al. 2007), and an initial temperature T c = 20 K. The average density, initial radius and 



free-fall time of the initial core is set by the initial volume average core surface density 

irrt 



(2) 



as 



9ttS 3 
16M 



M 



and 



ttM 



1/4 



(3) 
(4) 

(5) 



_64G 2 E 3 _ 

The initial core is placed in the center of a cubical simulation domain spanning 4 times the 
core radius i.e., (^domain = 4r c ) on each side, so that no part of the initial cloud except gas 
that is entrained into protostellar winds ever approaches the boundary. The initial core is 
immersed into a uniform ambient environment with density that is 0.01 times that at the 
edge of the initial core. Pressure balance between the core and environment is maintained 
by setting the temperature of the ambient gas to 100 times that at the edge of the initial 
core, T am b = 2000 K. The Planck mean opacity of the ambient gas is set to zero to ensure 
that it does not cool or radiatively heat the core. The cores are initialized with a turbulent 
velocity field chosen to put them in approximate balance between gravity and turbulent 
motions. Three Gaussian random fields are generated with power spectrum P(k) oc k~ 2 for 
the three velocity components, each normalized to have an integrated norm of unity over the 
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full spectral range sampled. We set the initial velocity in every cell equal to the components 
of the Gaussian random field times the one dimensional velocity dispersion, 



GM 



2(k p - l)r c 



G 2 MttS 



1/4 



(6) 



corresponding to the velocity at the surface of a singular polytropic sphere (McKee & Tan 



2003). The virial parameter a vir = ha 2 v GM/r c (Bertoldi fc McKee 1992) is 5 for k p = 3/2 



somewhat larger than the value of 15/4 in hydrostatic equilibrium ( McKee &: Tan|2003 ) . The 
kinetic energy is therefore initially larger than the gravitational energy, but a v - n decreases 
with time due to the decay of the turbulence. The initial radiation energy density is set to 
the value for a black-body radiation field with radiation temperature T r = 20 K everywhere. 

The earliest stages of high-mass star formation occur in infrared dark clouds (IRDCs) 



(Rathborne et al. 


2007 


) which are 


background ( 


Perault et al. 


1996; 


E 



of surface density in star forming regions from more tenuous sources with X r 
to more typical galactic star formation conditions with S ~ 1 g cm -2 , to E 



0.1 g cm -2 , 
^ 10 g cm -2 

(Beuther et al. 2002b Rathborne et al. 2006 Lopez-Sepulcre et al. 2010) or more in extra- 
galactic super- clusters ( |Turner et al.||200~0~)|McCrady fc Grahamp007l ) . Table |2~T| summarizes 



the parameters for each of the four computational models presented in this work. The initial 
conditions for these models have been chosen to study the collapse of galactic IRDCs with 
high but not atypical mass and surface density. Each of the initial simulation core states are 
rescaled versions of one another with identical density structure, virial ratio, velocity field 
and comparable peak resolution in every run. We have compared each of the simulations 
at equivalent time in units of the free-fall time of the initial cores. The homology between 
the runs is broken only by the presence or absence of outflows and by radiative effects. This 
choice of model parameters therefore probes the surface density dependence of radiative feed- 
back effects, and isolates the effects of protostellar outflows by holding all other parameters 
constant as they are turned on and off. 



We have largely followed the approach set forth by Krumholz et al. (2010) in choosing 
the parameters for the numerical simulations considered here. However, the simulations in 
this work differ from the earlier work in several ways. First, these simulations use an initial 
core mass of 300 M Q instead of 100 M . This choice was motivated by our desire to study the 
evolution of high-mass star systems and our expectation that protostellar outflows, which 
were not considered in the earlier models, will eject a significant fraction of the initial core. 
Secondly, the lowest initial surface density is E = 1.0 g cm -2 instead of E = 0.1 g cm -2 
as used in the earlier work. This choice is largely motivated by computational constraints. 
The high flow speeds present in runs with outflows necessitate smaller numerical time steps 
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than for non-outflow runs and thus increase the computational cost. A simulation with 
E = 0.1 g cm -2 would be particularly expensive due to its long free-fall time and the need to 
advance to these later times with numerical time steps that are limited by the cell crossing 
time of outflow-ejected gas. 

2.2. Refinement and Boundary Conditions 

The AMR capabilities of the code track the collapsing cores in three dimensions to grid 
scales of Axl ~ 25 AU on the finest AMR level. This resolution is achieved by discretizing 
the physical domain on the coarse onto a base grid of N$ cells. The placement of finer level 
grids up to the finest level L was determined by the refinement criteria that any gas denser 
than one half the density at the edge of the initial core be refined by at least one level. 
Further refinement is also triggered wherever the local Jeans number, J = ^jGpAx 2 / (7rcf) 
, exceeds 0.125, where Ax is the computational cell width on the coarser 
level, or wherever the local gradient of the radiation energy density |V 'E ra d\Ax / E ia <i exceeds 



Table 1. Simulation Parameters. 



E (g cm 2 ) 


1.0 


2.0 


2.0 


10.0 


wind 


on 


on 


off 


on 


M(M ) 


300 


300 


300 


300 


r c (pc) 


0.141 


0.100 


0.100 


0.0447 


fin (cm -3 ) 


7.3 x 10 5 


2.1 x 10 6 


2.1 x 10 6 


2.3 x 10 7 


a v /c s 


8.80 


10.5 


10.5 


15.6 


ts (kyr) 


50.7 


30.2 


30.2 


9.02 


-^domain (P*-0 


0.565 


0.400 


0.400 


0.179 


N 


128 


192 


192 


192 


max level 


5 


4 


4 


3 


Ax L (AU) 


28.4 


26.8 


26.8 


24.0 



Note. — Row 1: initial core surface density; Row 2: initial core mass; Row 3: initial core 
radius; Row 4: initial core velocity dispersion; Row 5: core free fall time; Row 6: linear size 
of the computational domain,; Row 7: number of cells per linear dimension on the coarsest 
level; Row 8: maximum refinement level; Row 9: computational resolution on the finest 
AMR level. 



(Trueloveet aL|1998 
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0.1. 

The simulations use a zero velocity gradient outflow boundary condition for the hydro- 
dynamics and Marshak boundary conditions for the radiation energy density with radiation 
flux at the edge of the computational domain for a 20 K background radiation field. The 
gravitational potential is specified at the edge of the domain as the sum of the multipole 
moments of the mass distribution in the computational volume as a function of time up to 
the quadrupole term. We adopt an equation of state with 7 = 5/3, appropriate for gas too 
cool for molecular hydrogen to be rotationally excited, but this choice is essentially irrelevant 
because the gas temperature is set almost purely by radiative effects. 



2.3. Optical Properties and Equation of State 

The radiation transport is handled by a frequency-integrated flux limited diffusion ap- 
proximation. We use the Planck and Rosseland mean dust opacities, Kp and respectively, 



of the Semenov et al. (2003) iron normal, composite aggregates model plotted in figure 
[TJ The simulations with protostellar winds introduce strong heating behind wind-driven 
shocks. When the thermal energy of the gas exceeds that of a gas with molecular weight 
0.6m p and temperature 10 4 K, we treat the gas as fully ionized with Rosseland opacity 
= 0.32 cm 2 g _1 , the value for Thompson scattering at solar metallicity. Our gray flux 
limited diffusion approximation cannot adequately represent the collisionally-excited line 
cooling processes that dominates at temperatures above the dust destruction temperature. 
However, it is critical to include them in order to ensure that shocked gas is able to cool. We 
therefore leave kj> = 10 -2 cm 2 g _1 for this gas, ensuring that it does not interact strongly 
with the ambient radiation field, but we also implement an approximate line cooling function 
A(T) to remove energy from gas above the dust destruction temperature and transfer it to 
the radiation field. We take A(T) from the function shown in figure 1 of Cunningham et al. 
(2006). In each time step, before we perform our ordinary flux limited diffusion radiation 
solve, we update the gas and radiation energy densities by implicitly solving the operator 
split-system 

^ = -(p/rfHT) (7) 
dE 

— = (p//i) 2 A(T) (8) 



using the LSODE ( Radhakrishnan & Hindmarsh 1993) Gear-type solver. In the above 



system, e is the specific thermal energy density of the gas, E is the radiation energy density, \x 
is the mean molecular weight and T is the gas temperature appropriate for a solar metallicity 
ionized gas mixture. 
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2.4. Protostellar Wind Model 



The ORION code includes a "star particle" algorithm to handle the formation of pro- 



tostars (Krumholz et al. 2004, 2007a). This algorithm provides for the creation of sub-grid 



star particles in those cells of the computation that become poised for gravitational collapse 
to spatial scales smaller than those that can be captured on the computational grid without 



spurious fragmentation (Truelove et al. 1997). The luminosity, radius and burning state 



of the star particle is advanced with the simulation according to the protostellar evolution 



model of McKee & Tan (2003) as updated by Offner et al. (2009). The protostellar evolu- 



tion model takes as input the mass and accretion history of the star as determined by the 
simulation, and as output predicts a protostar's radius and luminosity at any given time. 
The protostellar luminosity prescribed by this model enters the simulation as a source term 
in the radiation energy density equation, and the protostellar radius is used to compute the 
Keplerian velocity at the stellar surface, which affects outflows as described below. The 
protostellar luminosity prescribed by this model enters the simulation as a source term in 
the radiation energy density equation. 

The ORION star particle algorithm has been enhanced for this work to include the 
driving of bipolar outflows. Our outflow model is specified by the dimensionless parameters 
f w and f v , which set a wind launch speed as a fraction f v of the Kepler speed at the stellar 
surface and a mass flux that is a fraction f w of the rate of accretion onto the star or, 
equivalently, a fraction f w /(l + f w ) of the total mass that is either accreted onto the star 
or ejected in the wind. Since we are interested in the large scale impact of the protostellar 
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winds, we assume that the wind is injected over a range of radii determined by a function 
Xu)(l r l) with an angular dependence given by explicit forms for these functions are 

given below. The wind driving is imposed by operator split source terms in the gas density, 
momentum density and energy density equations with 

dp 



dt 
dpv 



dt 

dpe 

~dt 



where 



-f v v k ,i M u 



-M w4 x w (\*i\) Z(0i 



M(7 - 1) 



iGMi 



(9) 
(10) 

(11) 
(12) 



is the Keplerian speed at the surface of the star and r*^ is the protostellar radius; as remarked 
above, we use the value of r*^ given by the model of McKee & Tan (2003) as updated by 



Offner et al. (2009). For the simulations presented here we have set the wind-launched 
gas temperature as T w = 10 4 K, appropriate for an ionized wind. The corresponding rate 
of particle mass growth, particle wind mass ejection rate, acceleration, radial distance and 
spatial inclination of the i th star are 



Mi 



M„ 



l + /t 

UM 



-M. 



KKMOAi 



fu 



l + fu 



-M, 



KKMOAi 



9i 



VKKM04, 
X Xj, 

acos(ri • ji 



(13) 

(14) 

(15) 
(16) 
(17) 



where jj and Xj are the velocity and position of the i th particle, jj, Mkkmo4 and vkkmo^ are 
the sink particle angular momentum, accretion rate and acceleration for the case in which 



winds are absent as given by the algorithm of Krumholz et al. (2004) 



Values of the parameters f w and /„ can be both estimated from theory and constrained 



by observations. Theoretically, the X-wind (Shu et al. 1988) and disk wind (Pelletier & 
Pudritz 1992) models predict f w ~ 0.3, f v ~ 1 and f w ~ 0.1, /„ ~ 3, respectively. Both 
therefore suggest f w f v ~ 0.3. The total momentum P w carried by an observed outflow from 
a star of mass M* is related to and /„ by 

fwfv 



p 



(If 
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where Vk is the Keplerian velocity at the stellar surface. The peak of the stellar initial mass 
function is at M* « O.2M Chabrier (2005), and such stars typically have radii ~ 3R 
during their main accretion phases (e.g. model mC5H of Hosokawa et al. (2011)), so typical 
values of Vk are ~ 100 km s -1 . This value of Vk together with f w f v ~ 0.3 implies a net wind 
momentum flux of ~ 30 km s _1 per M of stars formed. 

A number of surveys have used measurements of P w and estimates of M* and v s , together 
with the relation given in equation (18), to constrain f w f v . Surveying the literature available 

0.3. More recent observational surveys of 

0.2 to 



as of 2000, |Richer et alT] ( |2000[ ) estimate /„,/, 

several nearby low-mass star forming regions indicate typical outflow momenta of 
~ 3.0 M km s _1 



Maury et al. 


2009 


Arce et al. 


2010 


Curtis et al. 


2010 



Ginsburg et al. 



2011). The physical properties of the driving sources of most of the surveyed outflows are 
not very well constrained. However, if we assume that the typical source has accreted half of 
its final mass M* ~ O.1M and has radius r* ~ 3R then equations (12) & (18) can be use 
to extract a range of outflow momentum parameters based on the results of these surveys of 
0.025 < f w f v < 0.38. 

Observationally, f w and f v can be better constrained from sources where observational 
measurements exist for both net outflow momentum and the net mass accreted onto the 
protostar M*. Curtis et al. (2010) have surveyed the outflow momentum in the young 



cluster NGC1333 and Sandell & Knee (2001) have estimated the mass of warm dusty gas in 
the collapsing envelope around the deeply embedded protostars that drive several of these 
outflows. In table 12.41 we list the intersection of those sources that have both dust mass 



measurements from Sandell & Knee (2001) and outflow momentum measurements by Curtis 



et al. (2010), excluding a few sources that Sandell & Knee (2001) indicated as near the 



edge of their field of view with unreliable flux densities, and excluding one tight binary 
source that could not be separately resolved in that work (IRAS 4A). We expect that the 
envelope masses of early class sources should be somewhat greater than the mass of the 
embedded protostar and we expect that the mass of later class I sources should exceed that 
of their envelope. In the absence of better mass constraints for the collection of class I 
and protostars in the present sample we adopt the assumption M* ~ M Dust as a rough 
approximation. We assume fiducial stellar radius of r* = 3 R following model mC5H of 



Hosokawa et al. (2011). From these assumptions, we can constrain f w f v from equations (18) 



& (12). The datum indicate a range of outflow launch parameters of 0.01 < f w f v < 0.15, 



with no strong statistical correlation of f w f v with the source spectral class. 

Wind launch speeds represent the Courant time-step constraint in a typical calculation, 
so large values of /„ impose a particularly onerous computational burden. We therefore 
choose a wind mass to stellar mass fraction on the high end of the theoretical guidance, 
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f w = 27% and the wind velocity parameter of, f v = 1/3. This yields a momentum flux 
injected by our wind model characterized by f w f v = 9%, which is toward the higher end of 
the observed range of rates of momentum injection by outflows from the low mass sources 
tabulated in table 12.41 

The function 

, . _ 1 T r~ 2 if 4Ax < r < 8Ax , . 

Xw{r) _-| Q otherwige 

is a normalized weighting kernel that determines the spatial scale of the wind injection where 
C\ is a normalization constant to the weighting kernel. The C\ is computed numerically in 
the ORION code so that numerical aliasing effects of the spherical wind injection region into 
the Cartesian grid are accounted for exactly. 

The remaining function £ describes the angular distribution of the wind mass and mo- 
Table 2. NGC1333 protostellar outflow data. 



Source 


Class 




p 


fwfv 


HRF42 





0.49 


0.058 


6.7xl0~ 4 


HRF43 


I 


0.36 


3.08 


0.057 


HRF44 





0.35 


3.17 


0.061 


HRF45 


I 


0.31 


0.28 


6.4xl0~ 3 


HRF46 





0.1 


0.44 


0.055 


HRF47 





0.24 


0.23 


7.8xl0~ 3 


HRF54 


I 


0.3 


0.10 


2.4xl0- 3 


HRF56 


I 


0.04 


0.11 


0.052 


HRF62 





0.32 


0.23 


0.0051 


HRF63 


I 


0.08 


0.07 


0.011 


HRF65 





0.07 


0.77 


0.17 


min 

mean 
max 








6.7xl0" 4 
0.039 
0.17 



Note. - Column 1: Hatchell et al. (2007) source number; Column 2: Spectral Class 



(Hatchell et al. 2007); Column 3: Progenitor mass (Sandell & Knee 2001); Column 4: Pro- 
genitor mass from the table B2 in the online supplementary data to (Curtis et al. 2010); 
Column 5: Implied outflow launch parameter assuming v w — 100 km s -1 . 
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mentum flux at the point where it is injected into the computational grid. We take this func- 



tion from Matzner & McKee (1999), who find that the momentum distribution of prestellar 



outflow injection asymptotically far away from the protostellar surface as a function of the 
polar angle 9 from the direction of the protostar's rotation is given by 



In 



(sm 2 9 + 9 2 n ) 



(20) 



where 8 is the so called "flattening parameter" that sets the opening angle of the wind. In 
the case of low-mass stars, Matzner & McKee (1999) suggest a fiducial value of 9 = 0.01 



and we adopt the same value here. While stars of type B or later direct momentum in a 



very well-collimated beam (Beuther et al. 2002a 2003 2004), O star winds are collimated 



somewhat more weakly, possibly due to the effect of ionization (Beuther & Shepherd 2005). 



Since we do not include ionizing radiation in these simulations, we do not attempt to model 
to model O star wind broadening. 



The large value of £ near 9 = in equation (20) requires particular care in implementing 



this model in a numerical code. We implement the driving function by averaging £ over the 
polar angle subtended by a grid cell A9 = atan(l/8) at the outer radius of the weighting 



kernel (equation (19)) as 



e+Ae/2 



£(Mc 




if | sin(| — 
otherwise. 



> 



Ax 



(21) 



Our choice to set £ to zero for angles close to ir/2 is driven by numerical considerations. If we 
allow the outflow to be injected into 4tt steradians around the star, its mass and momentum 
are sufficient to disrupt the early development of an equatorial disk. This behavior is an 
artifact of the necessarily poor resolution inside the wind launching region. We do not 
resolve the disk scale height, and this artificially puffs up the disk and reduces its mass and 
momentum density, rendering it far easier for the outflow to disrupt than it would be if its 
true scale height were resolved. We avoid this problem by reducing the outflow mass and 
momentum flux to zero in an equatorial belt that is at least one cell thick, ensuring that 
accreting material always has an uninterrupted path to the star. We note that Schonke~S; 



Tscharnuter (2011) have considered the effect of radiative and protostellar outflow feedback 



on the dynamics of the accretion disk in two dimensions at much higher resolution than 
the simulations considered here. Their simulations indicate that feedback effects can alter 
the accretion rate onto the star on shorter timescales and smaller length scales than have 
been resolved in this study. However, the emphasis of the present study is on the large scale 
radiative and feedback effects on the ambient core and we do not attempt to model this 
small-scale behavior. 
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The definite integral in equation (21) evaluates to 

9+A9/2 M 



-Ao/2 sin 2 # + #2 

1 



v ^2^T t an(^ + f ) , 
atan I — — I — atan 

#0 



ltan (9-f) 



On 



(22) 



The normalization constant 6*2(^0) — / $o)x«j( r )^ 3x is also computed numerically to 
exactly account for grid aliasing effects. Neglecting the grid aliasing effect, we find C 2 = 8.165 
by numerical integration. A second subtlety that arises in the numerical implementation 
of the wind driving is that care must be taken so that the momentum source terms impart 
exactly zero net momentum onto the star particles and the gas in the computational domain. 
If the position of the particle is allowed to vary continuously within its host cell the \ term 



in equation (19) may lead to an asymmetric driving. To overcome this problem, we bring 



the wind driving into symmetry with the numerical grid by rounding the particle position 
to the nearest half-integer multiple of the grid width 



Xi 2Ax nint 



25xJ 



(23) 



before computing the source terms (equations (13 10)) where nint is the "nearest integer" 
function. 

The wind momentum injection £ is significantly broadened at polar latitudes relative to 
the analytic prescription for £. Consequently less momentum is injected near the equator 
to satisfy the normalization constraint. To facilitate comparison of our numerical models 



with the analytic predictions of Matzner & McKee (2000) in £3.5 it is convenient to define 



the effective numerical flattening parameter O) eff over a range of angles that separate the 
outflowing gas from ambient gas: 

£(Mo = 0o,eff)=e(Mo = O.Ol). (24) 

We find that the above expression holds to with 10% for 9 > 10° with # ,eff = 5.75 x 10~ 4 . 



3. Results 

3.1. Large-Scale & Outflow Morphology 

The left column of figures [2] through [5] show the large scale evolution of each simulation 
from t = 0.2% tot = 0.8ig in increments of 0.2% . These plots show slices of density with the 
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color- mapping scaled by £ 3 / 2 following the scaling given in equation (pi). The spatial scale 
shown in each plot scaled by the initial core radius with each showing an area (2.5r c ) 2 . The 
slices are centered on the position of the primary protostar and oriented so that the angular 
momentum vector of the protostar is upwardly oriented on the page. As expected, once the 
scaling of cloud radius and surface density is accounted for, the regions that have not been 
penetrated by the outflow bow shocks collapse in a homologous manner with little dependence 
on the initial surface density. The protostellar outflows evacuate a shock-bounded cavity 
through the initial cores and the outflow cavities are the prominent features in the density 
slices in the left column by t ~ 0.3tff in the lower surface density cases and t ~ OAtg in 
the high surface density case. The propagation speed of the outflow bow shocks through 
the core and the width of the outflow cavities show a strong dependence on initial surface 
density. The lower surface density cores show significantly greater disruption due to the 
protostellar outflow feedback. This effect is due to 1) greater mechanical luminosity of the 
protostellar winds at lower surface density (an effect that we discuss in detail in £3.3 and 2) 
lower turbulent (cr 2 oc X 1 / 2 ) and thermal pressures (P c oc p c T c ~ £ 3 / 2 ) in the ambient cores 
which act to confine the propagation of the outflows. We note that the outflow-evacuated 
cavities in the cases with outflows emerge from the initial core toward the left side of the 
density slices in figure [5j The reason for this is that the primary star particle retains the 
momentum of the material that it accretes and therefore the star drifts away from the center 
of mass of the initial core with time (see £3.3). The outflowing material therefore emerges 
first from the thinnest side of the initial cloud, relative to the position of the primary star. 



3.2. Fragmentation & Star Formation 

The third column of figures [2] through [5] show the small scale evolution of the surface 
density. Each of the plots are scaled by the initial surface density E and the initial core 
radius, with each showing an area (0.1r c ) 2 centered on the position of the primary star. 
Deviations from the the homologous scaling with surface density exhibited on larger scales 
appear by t = 0.2t g. By t = QAtg in figure [3] and t = OMg in figure |4j notable differences in 
the disk around the primary protostar emerge. Progressing from the third row, showing the 
case with highest surface density, to the top row, showing the case of lowest surface density 
we note an increasing tendency of the disk around the primary star to fragment and generate 

. In the simulations with 
lower surface densities in the top two rows of figure |4j) we note a prominent increase in disk 
fragmentation in terms of the multiplicity of star particles and the presence of distinctly 
separate accretion disks around the primary and secondary protostars in comparison to the 
highest surface density case in the third row. 



spiral arms characteristic of a Toomre- unstable disk (Toomre||1964 
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The trend toward reduced disk stability and enhanced fragmentation with lower initial 



surface density is consistent with earlier analytic (Krumholz & McKee 2008) and numeri- 



cal (Krumholz et al. 2010) works that predict a threshold core surface density for sufficient 
radiative heating to inhibit disk fragmentation of £ ~ 1.0 g cm -2 , as well as with ob- 
servational data from infrared star-forming clouds that are consistent with this prediction 
( Lopez-Sepulcre et al.||2010 ). We will show in £3.4 that protostellar outflows provide a mech- 
anism for focusing radiative feedback in the poleward directions, away from the the infalling 



disk in the midplane, as predicted by Krumholz et al. (2005). Protostellar outflows should 



therefore raise the surface density threshold for high-mass star formation. The simulations 
presented here do not survey sufficiently low surface densities to quantify this effect, due 
to the computational cost of simulating protostellar winds inside low surface density cores 
with commensurately longer free fall times. Furthermore, we expect that the relative im- 
portance of this effect may depend on magnetic fields, which our simulations neglect, and 



their role in confining the outflow cavity (Hennebelle et al. 2011). We therefore defer more 



precise quantification of the effect of outflows on the minimum surface density threshold for 
high-mass star formation to future work that will include magnetic fields and survey lower 
surface density cores. 

We can isolate the effect of protostellar outflows on the small-scale evolution by com- 
paring the second and fourth rows of figures [2] through [5] which show the results of our 
numerical experiments with and without protostellar outflow ejection, both with the same 
initial surface density of £ = 2.0 g cm -2 . By t = OAtg enhanced radiation trapping in the 
case with outflows has lead to substantially warmer circumstellar gas in comparison to the 
case without outflows. The temperature structure in the £ = 2.0 g cm -2 case without proto- 
stellar outflows in the fourth row is more similar to the £ = 10.0 g cm -2 case with outflows 
in the third row than to the case at the same surface density with winds in the second row. 



As we will show in §3.4[ protostellar outflow cavities carve a path of reduced optical depth 
through the initial core that channel radiative flux away from the center of the core. This 
escape of radiative energy reduces the efficacy of radiative heating in the central regions of 
the collapsing core. 

Enhanced disk fragmentation associated with outflows and decreasing surface density is 
due to the reduced effectiveness of radiative heating. The right columns of figures [2] through 
[5] show the mass- weighted, column-averaged radiation temperature T r over the same regions 
as the column density projections in the center column. We note that the dense infalling 
gas is strongly coupled to the radiation field, so that the gas temperature T rj T r . Only 
the tenuous outflow-evacuated regions achieve sufficient post-shock temperatures to break 
the radiation-gas coupling by increasing the gas temperature beyond the dust destruction 
temperature. The column-averaged radiation temperature plots therefore allow us to probe 
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the temperature structure of the dense infalling gas without confusion from the projection 
of shock heated layers along the surfaces of the outflow cavity. The temperature structure of 
the gas on small scales shows even stronger dependence on the initial surface density than the 
column density. We note enhanced temperature with increasing surface density as early as 
t = 0.2tg when no stars are producing significant power via nuclear burning. The dependence 
of temperature on surface density at these times is solely due to the factors pointed out by 



Krumholz & McKee (2008); (1) Higher surface density cores have higher accretion rates, 
thereby generating higher accretion luminosity and (2) higher surface density cores have 
higher optical depth to more effectively trap radiation. The trend toward increasing gas 
temperature with increasing surface density is also present at all later times, due to additional 
radiative output from nuclear burning in the primary star in addition to enhanced accretion 
luminosity and enhanced radiative trapping. 

Stars with masses > 15M Q generate sufficient radiation pressure to exceed their grav- 
itational acceleration. The second column of figures [3j [4] and [5] show that by then the 
spherically-averaged radiation force from the central protostar exceeds the inward gravita- 
tion attraction acting on the dusty envelope of infalling gas in the case of S = 10.0 g cm -2 
and in the case of £ = 2.0 g cm -2 without outflows. In the case without outflows this strong 
radiative force drives the expansion of a bubble of circumstellar gas away from the central 
source. The early development of this bubble can be seen in the lower left panel of figure 
[3] at t = OAtff and the radial extent of the bubble grows to a size scale comparable to that 
of the initial core by t — 0.8tg as shown in figure |5j The radiation bubble emerges from the 
initial core on the left side of the density slices in figure [5j This is due to the drift of the 
primary star away from the center of mass of the initial core with time. The radiation bubble 
emerges first from the thinnest side of the initial cloud, relative to the position of the primary 
star. Accretion onto the primary star continues through the radiatively supported bubble 



via Rayleigh- Taylor unstable modes (Krumholz et al. 112009 Jacquet & Krumholz 2011) that 



develop dense, radiatively self-shielding spikes of infalling gas. The evolution of radiative 
bubbles in similar simulations without winds and without initial turbulence are discussed 



in detail by Krumholz et al. (2009). The sole difference between the radiation bubbles pre- 
sented here and those in the earlier work is that the bubbles presented here are considerably 
less symmetrical about the central source owing to the turbulent ambient environment. 

In the case with winds, the regions where the net force is dominated by the outward 
radiation force lie within the outflow cavity. These regions are dominated by outflow irre- 
spective of the radiation force. With protostellar outflow, in no case does radiation force 
exceed that of gravity acting on the infalling core gas. Consequently, no such radiation sup- 



ported bubbles form in any of the models with protostellar winds. As predicted in Krumholz 



et al. (2005), the cavities evacuated by protostellar outflows provide sufficient focusing of 
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the radiative flux in the poleward directions that accretion continues through the regions of 
the infalling envelope onto the disk that are not disrupted by the protostellar wind shocks, 
and the infalling motion of this gas is not interrupted by the effects of radiation pressure. 



3.3. Protostar Properties 

The upper left panel of figure [6] shows the time dependence of mass accretion onto star 
particles for each simulation and the middle left panel shows the time dependence of mass 
accretion onto the primary protostar. Abrupt jumps in the primary protostellar mass occur 
due to the merger of other particles that fall toward the primary star. The code has been 
constructed to merge star particles that cannot be resolved on the resolution scale of four 



computational zones on the finest level (see Krumholz et al. (2004)). Because these mergers 



are resolution-dependent effects, we have taken care to assure that the mass contribution to 
stellar sources due to mergers is small. As we have discussed earlier, the ambient core in 
the case of S = 1.0 g cm -2 is subject to the least efficient radiative heating and is the most 
susceptible to gravitational fragmentation and therefore the most difficult to resolve. The 
particle mergers account for ~ 15% of the total primary protostellar mass by the end of the 
simulation in the S = 1.0 g cm -2 case and about ~ 10% in the higher surface density cases. 
Most of the mass accumulated by merger events in the case of S = 1.0 g cm -2 occurs due to 
the merger of a secondary star particle of 3.OM at t — 0.52%. We do note, however, that 



Myers et al. (2011) have found that imposing a maximum mass threshold for mergers en- 
hanced fragmentation and limited the rate of mass growth of the primary star in more highly 
resolved but otherwise similar simulations. We therefore regard the stellar mass predictions 
in the models presented here as an upper limit. In comparing the cases with protostellar 
outflows we note that the total system mass in stars, M, exhibits a weak trend toward more 
rapid accretion with higher surface density even when the time is scaled by the free-fall 
time. In other words, Mts increases weakly with the surface density of the initial core, S. 
This is because outflows entrain and unbind less gas from the ambient core in cases with 
higher surface density. Therefore, the luminosity output from the more massive protostars 
is enhanced in the higher surface density cases. This contributes to more effective heating 
and decreased fragmentation of the ambient core in the higher surface density noted 



in £3.2 Consequently, the simulations at lower surface density fragment into low multiple 
systems earlier and are characterized by significantly reduced accretion to the primary star 
in comparison to the models at higher surface density. 

The upper right and middle right panels of figure [6] show the stellar luminosity and 
the speed of the protostellar wind driving from the primary protostar as a function of time. 
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Fig. 2. — Simulation graphics at t — 0.2tg for the parameters S = l.Og cm -2 , X = 2.0g cm -2 , 
X = 10. Og cm -2 , and £ = 2.0g cm -2 (without winds) from top to bottom. First column: 
p/E 3//2 on a (2.5r c ) 2 plane oriented so that the outflow launch direction lies in the plane of 
the image, pointing toward the top of the page. Black arrows indicate the velocity field. An 
arrow with length equal to 1/8 of the plot width indicates a flow speed of 100 km s -1 , and 
arrow lengths scale as vfv[. Second Column: Ratio of the radiation force magnitude to 
gravitational force magnitude. Third Column: Column density on a (0.1r c ) 2 plane aligned 
with the cardinal axes of the simulation, oriented so that the primary protostellar outflow 
direction is as close as possible to pointing vertically out of the page. Fourth Column: Mass- 
weighted radiation temperature projected in the same manner as the surface density in the 
third column. All plots are centered on the projected position of the primary star. White 
markers indicate star particles. 
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Fig. 3. — Same as figure [2] but at t = 0.4%. 
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4. — Same as figure [2] but at t = 0.6% 
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Fig. 5. — Same as figure [2] but at t — 0.8tg. Only the surface density parameter cases of 
£ = 2.0g cm -2 , X = 10. Og cm -2 , and S = 2.0g cm~ 2 (without winds) were run to this time. 



Nuclear burning dominates the radiative output from stars with mass exceeding 5M Q and 
the primary protostar dominates the stellar feedback into the core. The rapid increase in 
mechanical feedback from winds in all of the models is driven by the Kelvin-Helmholtz 



contraction of protostar (Shu et al. 1987). Stellar contraction causes the escape speed at 



the protostellar surface to increase which in turn increases the wind ejection speed with 
: fesc/3 (see £2.4). By t = 0.2tg all of the models have undergone sufficient contraction 



v 



to ignite deuterium burning in the stellar cores, leading to a rapid increase in luminosity. 
Rapid contraction of the stellar surface continues as deuterium is exhausted in the prestellar 
core, giving rise to a commensurate increase in wind speed from t = 0.2tg to t — 0.3tg. 



Contrary to the increase in the effects of radiative feedback with surface density, we 
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find that the effects of mechanical feedback on the ambient cores decrease with the initial 
surface density. For t > 0.4tg the primary wind speeds show a noticeable decrease with 
increasing surface density. This occurs due to a mismatch between the Kelvin time for 
the contraction of the stellar surface and the free-fall time of the ambient molecular cloud 
core. The former depends mostly on the protostellar mass whereas the latter scales as 
t$ oc £ -3//4 (equation (JH])). This means that lower surface density simulations form stars 
that undergo more rapid Kelvin contraction per unit free-fall time and thereby eject more 
powerful winds at equivalent stages of collapse. We note that this result is driven by the 
assumption built into our numerical model that the wind ejection speed is proportional 
to the escape speed at the protostellar surface. As discussed in §3.2 , the more powerful 



winds contribute to the enhanced disruption of the ambient core in the lower surface density 
simulations. However, the overarching trend toward enhanced disruption of the ambient core 
with lower surface density would remain even if the wind speed as a function of core free-fall 
time were independent of surface density. First, the outflow break-out time from the core, 
^outflow = r c/v w would scale, relative to the free-fall time, as £ outflow As ~ S 1//4 , resulting in 
later outflow emergence from higher surface density cores. Second, higher surface density 
cores enhance the confinement of outflow cavities due to higher ambient pressures. 

The cascade of turbulent motions through the collapsing cloud strongly affects the mo- 
tion of the primary star as a function of time. The lower left panel of figure [6] plots the 
distance between the primary star, at position r p , and the center of mass of the system, at 
position ycom- The drift of the primary star away from the center of the system is a result 
of the random velocity field. Because most of the turbulent energy is in large wavelength 
modes, the denser center of the cloud will tend to have a different velocity than the lower 
density edges. The star initially forms in a compressional mode that is much larger than the 
local reservoir of gas that starts the protostar, but smaller than the diameter of the core. 
This initial reservoir of gas that starts the star is not co-moving with the center of mass 
of the core. Therefore, while the star does initially form at the bottom of the gravitational 
potential well it has a finite kinetic energy and can oscillate within the well. As a result the 
location of the star drifts from where it first forms in our simulations. 

The cascade of turbulent motions through the collapsing cloud also strongly affect the 
orientation of the primary star as a function of time. The lower right panel of figure [6] 
shows the angle of inclination of the angular momentum accreted by the primary protostar 
particle relative to the z direction. We note that the total angular momentum of the initial 
cloud has an orientation that is inclined 27° to the z direction. Because no radiative heating 
feedback is present in the initial state of the cloud, the early evolution of the simulations 
are characterized by fragmentation of the densest portion of the initial cloud into 2 to 4 
gravitationally bound particles, as shown in figure [2j By t = 0.3ts, the radiative feedback 
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from the primary star is sufficient to suppress local fragmentation and the initial fragments 
merge. The angular momentum of the primary particle at early time varies over ~ 90° 
as a consequence of variation in the angular momentum of the surrounding gas and as a 
consequence of the coalescence of the initial fragments. After t = 0.3ig the primary star has 
built up a sufficient moment of inertia relative to the rate of angular momentum deposition 
by accretion that the orientation of the star changes less rapidly. Subsequent evolution 
(t = 0.3tff tot = 0.8% ) of the orientation of the primary star is characterized by a rotation of 
20° to 40°. We note that the change in angular momentum in the lower surface density cases 
(X = 1.0 g cm -2 and £ = 2.0 g cm -2 ) occurs in abrupt jumps, whereas the higher surface 
density case is less susceptible to numerical fragmentation and therefore are characterized 
by relatively smoother change. 

The cumulative distribution of stellar masses at t = 0.5/% for each of the runs shown in 
figure [7j Both the highest surface density case with X = 10.0 g cm -2 and the case without 
outflows at surface density S = 2 g cm -2 collapse to a single star of 70% of the total mass 
accreted, consistent with the qualitative similarity in the small-scale temperature structure 



noted in £3.2 The lower surface density cases on the other hand fragment into binary 
systems with > 1 M Q secondaries present by t = OAts- These results demonstrate that the 
absence of outflows and/or higher initial surface densities result in less fragmentation and 
the production of fewer, more massive stars. 



3.4. Radiation Focusing 

The radiative feedback from stellar sources into the ambient cloud is not isotropic. Four 
possible causes of anisotropic stellar output include (1) the clearing of an optically thin 
outflow-swept cavity along the poles of the star, (2) shielding in the midplane due to the 
presence of an optically thick accretion disk, (3) non-uniform optical thickness of the ambient 
cloud due to the displacement of the star relative to the center of mass of the ambient cloud, 
and (4) non-uniform optical thickness of the ambient cloud due to the turbulent density 
structure. In this section we consider the distribution of radiative output from the primary 
star in each simulation as a function of solid angle to elucidate the importance of each of 
these effects in shaping the radiative feedback into the ambient cloud. 

Figure [8] shows the radiative flux from the primary protostar as a function of spherical 
angle in each model at times t = 0.3, 0.4, and 0.6%, normalized to the radiative flux that 
would be expected in an isotropic environment F isotropic • r = L p j (47rr 2 ). The plots shown in 
figure [8]have been constructed from slices of the radial radiation flux at a distance of 1500 AU 
from the primary star, rotated into the coordinate system where the angular momentum of 
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Fig. 6. — Stellar properties as a function of time for each set of simulation parameters. Upper 
Left: Total mass in stars. Upper Right: Primary star luminosity. The luminosity has been 
smoothed using a 200 yr moving average to eliminate the high frequency contribution of the 
accretion luminosity in this plot. Middle Left: Primary star mass. Middle Right: Primary 
protostellar wind speed. Lower Left: Position of the primary star relative to the center of 
mass of cloud. Lower Right: Angle between the primary star's angular momentum vector 
and the z-axis. 



- 25 - 



1.0 



1.0 



2.0 



0.8 



E = 10.0 

E = 2.0, no wind 



C-0.6 

V 



0.4 



0.2 



0,0, 




10° 



10 : 



M[M G ] 



Fig. 7. — Fraction /(< M) of total stellar mass contained in stars with mass < M as a 
function of the total mass in stars for each of the runs at time t = 0.5tg. 

the primary star points upward. The angular distribution of the radiative flux is relatively 
insensitive to the position of the spherical slice, provided that the slice has radius greater 
than that of the disk around the primary protostar. In varying the radius of the spherical 
slice from 1500 AU to 10 4 AU, we note a decreased contrast between regions of low and 
high radial flux by about ~ 25%, which we attribute, in part, to the flux limited diffusion 
approximation used for the radiation transport. However, the location and extent in solid 
angle of the large-scale features is independent of the position of the slice. The azimuthal 
coordinate facing the thinnest edge of the cloud due to the motion of the star relative to 
the center of mass of the system is centered in each of the plots. Regions of peak outward 
radiative flux correlate very well with outflow ejection, and regions of low outward radiative 
flux correlate very well with the presence of the accretion disk near the midplane. The 
clearing of low optical depth paths of escape by outflow ejection and shielding by the dense 
accretion disk in the midplane are therefore the dominant effects in focusing the radiative 
feedback from the star. As gas falls in toward the primary protostar, it carries its angular 
momentum with it. In the highest surface density gas the accretion flow is fairly smooth 
because radiative heating raises the pressure near the primary star and prevents gas around 
it from clumping up. At lower surface densities, however, the accreting gas may be partially 
collapsed under its own gravity, or may even have collapsed completely to form stars that 
then merge with the primary. As a result, angular momentum tends to be accreted in 
distinct lumps, leading to rapid reorientation of the accreting star over short timescales. 
The radiative flux is far more focused as bipolar in the case of £ = 10.0 g cm" 1 , consistent 
with the narrower geometry of the outflow cavity in this case. 



-26- 



Figure [9] shows the azimuthally-averaged radiative flux from the primary protostar as 
a function of polar angle in each model at times t = 0.3, 0.4, and 0.6tfj. At each of these 
times, the effect of the presence of protostellar outflow cavities is clearly evident showing polar 
radiation flux 1.7 to 15 times that at the midplane. The models with protostellar outflow 
show enhancement of the polar flux relative to the case without outflows by comparable 
factors. The degree of poleward focusing of the radiation flux diminishes with time due to 
broadening of the outflow evacuated cavities that focus the radiation. We note that polar 
flux in the £ = 1.0 g cm -2 and X = 2.0 g cm -2 cases at t — OAtg underrepresents the flux 
focusing due to outflow cavities because the outflow cavities are tilted by ~ 20° relative 
to the orientation of the primary protostar shown in figure [8j Similar radiation focusing 
effects were also shown in Krumholz et al. (2005), where the authors considered the effects 
of the presence of outflow cavities on radiation escape from the infalling envelope around 
massive protostars. Using static radiative transfer models they showed that focusing of the 
radiative force from the central star throughout the outflow cavity results in a reduction of 
the equatorial radiative flux relative to a control model without outflows by factors of 1.7 
to 14, depending on the width and shape of the region of low optical depth in the outflow 
cavity. Therefore, the models presented here support the Krumholz et al. (2005) prediction 
that outflows reduce the Eddington radiation pressure barrier to high-mass star formation 
by reducing the radiation force exerted in the infalling cloud gas. However, it has also been 
shown that fully 3D Rayleigh- Taylor modes can remove the Eddington barrier even when 
protostellar outflows are neglected (Krumholz et al. 2009). 



3.5. Star Formation Efficiency 



Matzner & McKee (2000) give an analytic model for the core to star formation efficiency 

e core = 1 - MJM, (25) 

where M e j is the net mass ejected from the core by entrainment into the protostellar outflow 
and M is the initial core mass. For the case of an unmagnetized core the model predicts 

2X(c„) , s 

(26) 



1 + ^1 + 4(1 + ^)2X2^)' 

where the dimensionless function X is defined by 

x = c ^(f)tt^ (27) 

\ y 0/ JwV w 

v csc is the escape speed from the edge of the core, and the parameter c g depends on the 
core density profile and free-fall time relative to the effective timescale for the wind driving. 
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Fig. 8. — Radial component of the radiation flux, normalized to the isotropic flux at 1500 AU 

from the primary protostar. Columns indicate the times t = 0.3tg, t = OAtg and t = 0.6tg 

from left to right and the rows indicate the simulation parameters of £ = 1.0 g cm~ 2 , 

£ = 2.0 g cm -2 , S = 10.0 g cm -2 and £ = 2.0 g cm -2 without winds from top to bottom. 

The coordinate system is defined such that the angular momentum of the primary star points 

northward and the azimuthal coordinate facing the thinnest edge of the cloud due to the 

motion of the star relative to the center of mass of the system is centered in each of the plots, 
th 

Contours of the 75 percentile column density from r = to 1500 AU are shown as dashed 
lines and contours of the radial velocities of 20 km s -1 and 50 km s _1 at r = 1500 AU are 
shown as solid lines. 
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Fig. 9. — Radial component of the radiation flux at 1500 AU in radius from the primary 
star, normalized to the isotropic flux, as a function of polar angle. The flux shown at a given 
9 is a volume average over a pair of rings at polar angles 9 = and 180° — 9 that cover 
all azimuthal angles <fi. The coordinate system is oriented so that 9 = corresponds to the 
rotation axis of the primary star and the direction in which the wind is launched. 
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For an unmagnetized core with profile p oc r kp , equation (A19) of Matzner & McKee (2000) 
provides an estimate of 

tt(9-3/c p )(4-A; p ) ft, 



c q>>1 = v p, \ p> - (28) 

for steady winds. The estimate depends on the the age of the steady wind, t w and is valid 
for c g ^> 1. In the limit of an impulsively driven wind, t w — > 0, Matzner & McKee (2000) 
provide an estimate of 

c = h^h. (29) 

To compare our simulations with the analytic model, we choose c g by interpolating between 
the two limits as: 

_ 7r(9-3fc„)(4-fc„) ft w \ / 9-3fe p 
C9 " 8^ UJ + V 
For fc p = 3/2 this expression becomes 



c g = 2.52— + 1.13. (31) 

The mass-weighted average wind speed that characterizes the wind momentum injection into 
the core is 

v w = ^2 TU I M w ,iV wA dt. (32) 

stars JwMi J 



Values of v w and v esc are given in table 3.5 Because our numerical wind injection approach 



is based on volume-averaged quantities inside of an 8-cell radius wind injection sphere as 



described in £2.4, the effective numerical flattening parameter is 6o, e s — 5.75 x 10 -4 for winds 
with an opening angle > 32° and we shall use this as the flattening parameter for the purposes 
of comparing the simulations to the analytic model. We note that X depends logarithmically 



on the flattening parameter (equation (26)) and therefore the model prediction is not very 



sensitive to the estimate of the effective numerical flattening angle. 

Due to constraints on computational time, we have not run numerical simulations suffi- 
ciently long to determine the final star formation efficiency. To facilitate comparison between 
the numerical and analytic model, we focus our attention to the ratio of the mass ejected by 
winds to the total stellar mass, 

_ j^gj _ 1 ~ e corc /„„x 
fcwind — ^ M ~ c ' ^ ' 

Z^stars 1Vli e core 

where M e j is the total mass ejected from the system. This quantity can be computed as a 
function of time throughout the simulation. For the purpose of comparing our numerical 
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Fig. 10. — Left: Mass with outward radial speed greater than the escape speed of the system 
as a function of the total mass in stars. Right: Mass with outward radial speed greater than 
the escape speed relative to the total mass in stars as a function of time. 

simulations to the analytic prediction, we heuristically define the ejected mass as any mass 
that has been either ejected from the simulation domain or that is propagating with a 
sufficient radial component of velocity away from the center of mass of the system to overcome 



its gravitational binding to the system. The left panel of figure [10] shows a plot of the total 
wind mass ejected in each simulation as a function of the total mass in stars and the right 
panel shows the simulation result for the wind ejection efficiency e w i n d = M e j/ J^stars ^ as a 
function of time for each simulation. We note that wind-launched gas in the highest surface 
density simulation with S = 10.0 g cm -2 emerges from the initial core relatively late in the 
simulation at t ~ 0.6%. At the end of the simulation, much of the wind-launched gas is still 
entrained in the limbs of the outflow cavity. This is a transient effect that is not included 
in the analytic model and we therefore will focus our attention in comparing the analytic 
prediction to only the S = 1.0 g cm -2 and £ = 2.0 g cm -2 models. By inspection of figure 



10 , we adopt the values of t w = 0.5% and t w = 0.3% as characteristic of the age of the winds 



in the £ = 1.0 g cm 2 and £ = 2.0 g cm 2 simulations respectively. 



The analytic predictions given by equations (26) & (33) for the outflow mass- weighted 
wind speed v w for each simulation are given in table 3.5 The wind ejection efficiency at the 
end of each simulation (t = t end ) are listed in the table as e win d, simulation for comparison to 
the model predictions. 



In drawing conclusions from comparing the analytic model result for the wind ejection 
ratio e win( j to the simulation result e win d, simulation it is important to bear in mind that the 
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analytic model predicts the fraction of the initial core that is ejected by a wind over the 
entire course of evolution of the initial core, whereas the numerical model is only averaged 
over the duration of the core evolution up to the end of the simulation. Furthermore, the 
analytic model cannot account for the time dependence associated with the propagation of 
protostellar wind shocks through the highly inhomogeneous ambient core. 

For the lowest surface density case £ = 1.0 g cm" 2 , the simulation has ejected more mass 
with the wind per unit of mass accretion than predicted by the analytic model by a factor of 
2.3. This is likely due to efficient entrainment of core gas into the wind via the interaction 
of the protostellar winds as they propagate through the network of dense filaments in the 
ambient core at early time. By t = 0.4 tg, a wide wind cavity has been cleared through 
the initial core in the simulation (see the top row of figure [3]). The time-integrated outflow 
ejection behavior in the simulation is biased toward this early-time behavior because the 
simulations are only advanced to t = 0.6tg in the E = 1.0 g cm -2 case and t = 0.8tg in the 
other models. In figure 10, we note that evolution of the system past 0.5tff carries forward 
with less efficient entrainment of core gas into the outflows as the outflowing gas at later time 
propagates unimpeded through the wind channel. We expect that with continued evolution 
the system would asymptote to a steady state value of e cor c that is closer to the model 
prediction. However, constraints of computational cost have prevented us from testing this 
expectation. 

The intermediate surface density case £ = 2.0 g cm -2 exhibits low ejection efficiency 



Table 3. Outflow ejection. 



£ (g cm -2 ) 


1.0 


2.0 


10.0 


^end (tfi) 


0.6 


0.8 


0.8 


v csc (km s" 1 ) 


4.27 


5.08 


7.60 


v w \t=t end (km s" 1 ) 


87.7 


72.2 


71.0 


^core 


0.70 


0.73 




Cwind, simulation 


1.06 


0.342 


0.0563 


Cwind 


0.42 


0.370 





Note. — Simulation results (rows 1-5) and analytic model predictions (rows 6 & 7). The 
columns indicate the cases of £ = 1.0 g cm -2 , £ = 2.0 g cm -2 and E = 10.0 g cm -2 from 
left to right. As discussed in the text, the E = 10.0 g cm" 2 simulation was not evolved 
sufficiently far in time to compare with the analytic model. 
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due to during the initial core-crossing of the outflow from t — to t — 0.3ig. The ejection 
efficiency plateaus at later time with e core varying between 0.34 and 0.42 at later time. We 
note that by the end of the simulation at t = 0.8tg, the outflow ejection efficiency in the 
simulation is 92% of the analytic prediction. 

We note that the star formation rate per free-fall time in our simulations, eg = e cove tff/t en d, 
is much greater than the observationally-estimated value of a few percent found by |Krumholz| 



&Tan 


(2007 


) and 


Evans et al. 


(2009) 



both observationally and theoretically. 



On the observational side, the Krumholz & Tan (2007) and Evans et al. (2009) estimates 



were for gas clumps at densities of at most ~ 10 5 cm -3 , and the typical objects at this density 
are ~ 10 4 M parsec-sized clumps that are forming entire star clusters. In comparison, the 
presetellar cores that we have simulated are much smaller and denser: n ~ 10 6 cm" 3 , 
r ~ 0.1 pc, M ~ 10 2 Mq. In the terminology of McKee & Ostriker (2007), they are "cores" 



rather than "clumps" . There are no observational measurements for the value of eg in such 
structures. Indeed, Krumholz & Tan (2007) commented that eg must reach values ~ 1 
rather than ~ 0.01 at some density higher than what then current observations probed. 
Furthermore, it is clear that 100 M & cores forming massive stars must be rare exceptions, 
since massive stars are rare. Since measurements of eg only provide statistical averages, it is 
possible that a few n ~ 10 6 cm -3 cores like the ones we have studied undergo rapid collapse, 
but that there are more numerous structures at similar density that are not undergoing rapid 
monolithic collapse, so that the average value of eg is much lower than in the core we have 
simulated. 

On the theoretical side, values of eg ~ 0.01 are expected only in regions where there 



is turbulence at roughly virial levels (Krumholz et al. 2005). In our simulations, while we 



start out with such turbulence, this decays rapidly. Since these simulations contain a single 
massive star with a dominant outflow, there is nothing to drive turbulence in the core, and 
this allows eg to rise rapidly as the turbulence becomes sub-virial. One can see this effect 
in figure 6: the total mass in stars rises very slowly at first, and accelerates as time goes on 
and the turbulence decays. 



4. Summary 

We report the results of several AMR radiation-hydrodynamic simulations of the collapse 
of massive star forming clouds using the ORION code. These simulations are the first to 
include the feedback effects of protostellar outflows, radiative heating, and radiation pressure 
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in a single computation. In these simulations, the initial density profile, velocity spectrum, 
virial ratio and numerical resolution are held constant. The simulations are scaled to different 
surface density to study the environmental dependence of the outflow and radiation feedback, 
and in one case the surface density is held constant but outflow feedback is turned off to 
isolate the effect of protostellar outflow. 

Comparison of models with protostellar outflow feedback and surface densities S = 
1.0 g cm -2 , £ = 2.0 g cm -2 and £ = 10.0 g cm -2 at equivalent free-fall time shows that the 
higher surface density clouds exhibit enhanced radiative heating feedback, diminished disk 
fragmentation and host more massive primary stars with less massive companions. However, 
the effects of outflow feedback diminish with increased surface density. Lower surface den- 
sity clouds have longer free-fall time and therefore undergo more Kelvin contraction in the 
primary protostellar core, leading to more powerful outflows and more effective mechanical 
feedback. Furthermore, lower surface density clouds give rise to protostellar outflows with 
shorter core crossing time relative to the core free-fall time and these clouds are consequently 
influenced by the effect of outflows at relatively earlier stages of collapse. 

Comparison of models with and without outflow feedback at surface density X = 
2.0 g cm -2 indicates a strong coupling between outflow and radiative feedback on the parent 
cloud. Outflow activity produces polar cavity of reduced optical depth through the ambi- 
ent core. Radiation focusing in the direction of outflow cavities is sufficient to prevent the 
formation of radiation pressure-supported circumstellar gas bubbles, in contrast to models 
which neglect protostellar outflow feedback. With outflows, the radiative flux in the pole- 
ward direction is enhanced by 1.7 to 15 times that in the midplane. Sheets with outward 
radiative flux reduction up to an order of magnitude appear near the equatorial latitude 
of the primary star in all of the models with protostellar outflow. This result is consistent 



with the predictions of Krumholz et al. (2005) that focusing of the radiative flux from the 



central star throughout the outflow cavity results in a reduction of the equatorial radiative 
flux relative to a control model without outflows by factors of 1.7 to 14, depending on the 
geometry of the outflow cavity. As a result the radiative heating and outward radiation force 
exerted on the protostellar disk and infalling cloud gas in the equatorial direction is greatly 



diminished by the presence of the outflow cavity, and our models support the Krumholz 



et al. (2005) prediction that outflows reduce the Eddington radiation pressure barrier to 
high-mass star formation by reducing the radiation force exerted in the infalling cloud gas. 
Precisely determining the effect of outflows on the threshold density prediction for massive 
star formation will require examination of simulations at lower cloud surface density than 
the models presented here. Furthermore, we expect that the relative importance of this 
effect may depend on the role of magnetic fields in confining the outflow cavity. Future 
works should therefore examine the effect of protostellar outflows in lower surface density, 



magnetized clouds. 
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